clear all;
close all;
clc;
%% PLEASE MODIFY THIS PART
% W={'0.1','0.2','0.4','0.6','0.8'};
% w=[0.1:0.1:2.0];
w = [0.1:0.1:2.0];
Mode_and_LatticeSize = [20001,40000];
% W={'10'};

%% SHOULD NOT BE NECESSARY TO CHANGE
%% FITTING STAGE 1
for i=1:length(w)
    % Load data ============================
    % w(i)=str2num(W{i});
    a=load(['C:\Documents and Settings\Xinyu\Desktop\LPCVD\W_control\Simulations\Open_loop\KMC_1\W_' num2str(w(i)) '_A_0.0\Output_r2_1.dat']);
    b=load(['C:\Documents and Settings\Xinyu\Desktop\LPCVD\W_control\Simulations\Open_loop\KMC_1\W_' num2str(w(i)) '_A_0.0\Output_m2_1.dat']);
    time(i,:)=a(:,1)';
    rough2(i,:)=a(:,2)';
    slope2(i,:)=b(:,2)';

    % Fitting ==============================
    % time(i,:)=load(['1d_en0_T480_W' W{i} '_L4e4_t5e2_t.txt']);
    % rough2(i,:)=load(['1d_en0_T480_W' W{i} '_L4e4_t5e2_r2_1.txt']);
    % slope2(i,:)=load(['1d_en0_T480_W' W{i} '_L4e4_t5e2_m2_1.txt']);
%    cost_fun=@(para) sum(0*(model_r2(para,time(i,:))-rough2(i,:)).^2+1*(model_m2(para,time(i,:))-slope2(i,:)).^2);
    cost_fun=@(para) sum(0*(model_r2_mean_1D(para,time(i,:),Mode_and_LatticeSize)-rough2(i,:)).^2+1*(model_m2_mean_1D(para,time(i,:),Mode_and_LatticeSize)-slope2(i,:)).^2);
    para=fmincon(cost_fun, [1000 .01],[0 0],[0]);
    p(i,:)=para;
    r2(i,:)=model_r2_mean_1D(p(i,:),time(i,:),Mode_and_LatticeSize);
    m2(i,:)=model_m2_mean_1D(p(i,:),time(i,:),Mode_and_LatticeSize);
    
    % Verification =========================
    figure(1);
    plot(time(i,:) , rough2(i,:), time(i,:), r2(i,:),'--');
    figure(2);
    plot(time(i,:) , slope2(i,:), time(i,:), m2(i,:),'--');
    % plot(time(i,:) , rough2(i,:), time(i,:), r2(i,:),'--');
    hold on;
end

%% FITTING STAGE 2
% h=figure;
% plot(time , rough2*0.04, time, r2,'--');
% h=figure;
% plot(time, slope2, time, m2,'--');
% 
% h=figure;
% plot(w',p(:,1),'rs');
% % fit_c2=polyfit(w',log(p(:,1)),4);
% fit_c2=polyfit(w',p(:,1),4);
% hold on;
% % plot(w',exp(fit_c2(1)*w.^4'+fit_c2(2)*w.^3'+fit_c2(3)*w.^2'+fit_c2(4)*w.^1'+fit_c2(5)),'--r');
% plot(w',fit_c2(1)*w.^4'+fit_c2(2)*w.^3'+fit_c2(3)*w.^2'+fit_c2(4)*w.^1'+fit_c2(5),'--r');
% hold off;
% 
% h=figure;
% plot(w',p(:,2).*p(:,1),'b^');
% fit_sigma2=polyfit(w',p(:,2).*p(:,1),4);
% hold on;
% plot(w',fit_sigma2(1)*w.^4'+fit_sigma2(2)*w.^3'+fit_sigma2(3)*w.^2'+fit_sigma2(4)*w.^1'+fit_sigma2(5),'--b');
% hold off;